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Abstract 

A  general  formulation  for  a  comprehensive  fuel  cell  model,  based  on  the  conservation  principle  is  presented.  The  model  formulation 
includes  the  electro-chemical  reactions,  proton  migration,  and  the  mass  transport  of  the  gaseous  reactants  and  liquid  water.  Additionally,  the 
model  formulation  can  be  applied  to  all  regions  of  the  PEM  fuel  cell:  the  bipolar  plates,  gas  flow  channels,  electrode  backing,  catalyst,  and 
polymer  electrolyte  layers. 

The  model  considers  the  PEM  fuel  cell  to  be  composed  of  three  phases:  reactant  gas,  liquid  water,  and  solid.  These  three  phases  can  co-exist 
within  the  gas  flow  channels,  electrode  backing,  catalyst,  and  polymer  electrolyte  layers.  The  conservation  of  mass,  momentum,  species,  and 
energy  are  applied  to  each  phase,  with  the  technique  of  volume  averaging  being  used  to  incorporate  the  interactions  between  the  phases  as 
interfacial  source  terms.  In  order  to  avoid  problems  arising  from  phase  discontinuities,  the  gas  and  liquid  phases  are  considered  as  a  mixture. 
The  momentum  interactions  between  the  fluid  and  solid  phases  are  modeled  by  the  Darcy-Forchheimer  term.  The  electro-oxidation  of  H2 
and  CO,  the  reduction  of  02,  and  the  heterogeneous  oxidation  of  H2  and  CO  are  considered  in  the  catalyst  layers.  Due  to  the  small  pore  size 
of  the  polymer  electrolyte  layer,  the  generalized  Stefan-Maxwell  equations,  with  the  polymer  considered  as  a  diffusing  species,  are  used  to 
describe  species  transport. 

One  consequence  of  considering  the  gas  and  liquid  phases  as  a  mixture  is  that  expressions  for  the  velocity  of  the  individual  phases  relative 
to  the  mixture  must  be  developed.  In  the  gas  flow  channels,  the  flow  is  assumed  homogeneous,  while  the  Darcy  and  Schlogl  equations  are  used 
to  describe  liquid  water  transport  in  the  electrode  backing  and  polymer  electrolyte  layers.  Thus,  two  sets  of  equations,  one  for  the  mixture  and 
another  for  the  solid  phase,  can  be  developed  to  describe  the  processes  occurring  within  a  PEM  fuel  cell.  These  equations  are  in  a  conservative 
form,  and  can  be  solved  using  computational  fluid  dynamic  techniques. 

©  2004  Elsevier  B.Y.  All  rights  reserved. 

Keywords:  PEM  fuel  cell;  Mathematical  modeling;  Two-phase  flow 


1.  Introduction 

The  use  of  fossil  fuel  is  a  major  source  of  air  pollu¬ 
tion  and  contributes  to  global  warming.  The  transportation 
sector  is  a  major  consumer  of  fossil  fuel,  thus  eliminating 
or  reducing  pollution  from  transportation  sources  is  a  ma¬ 
jor  policy  objective.  Polymer  electrolyte  membrane  (PEM) 
fuel  cells  convert  the  chemical  energy  of  H2  and  O2  di¬ 
rectly  into  electrical  energy,  with  water  and  heat  as  the 
only  by-products.  The  low  operating  temperature  allows  for 
quick  start-up,  and  the  high  power  density  and  mechanically 
robust  construction  make  PEM  fuel  cells  an  attractive  re- 
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placement  for  the  internal  combustion  (I.C.)  engine.  Cur¬ 
rently,  PEM  fuel  cells  are  not  a  commercially  viable  alter¬ 
native  to  the  I.C.  engine;  however,  a  greater  understanding 
of  the  processes  occurring  within  a  PEM  fuel  cell  can  aid  in 
commercialization  by  providing  power  output  and  efficiency 
gains. 

Insight  on  PEM  fuel  cells  can  be  obtained  through  math¬ 
ematical  modeling.  Mathematical  models  use  fundamental 
equations  to  simulate  the  processes  occurring  within  a  PEM 
fuel  cell.  Although  several  processes  occur  within  a  PEM 
fuel  cell,  three  key  processes  have  the  greatest  impact  on 
PEM  fuel  cell  performance:  the  electro-chemical  reactions 
in  the  catalyst  layers,  proton  migration  in  the  polymer  elec¬ 
trolyte  membrane  layer,  and  mass  transport  within  all  regions 
of  the  PEM  fuel  cell.  These  processes  have  been  addressed  by 
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Nomenclature 

a 

activity 

A 

area  (m2) 

An 

reactive  area  per  volume  (m-1) 

b 

body  force  (N  m-3) 

b f 

ratio  of  forward  to  backward  reaction  rate  con¬ 
stant  (mole  m-3  or  m3  mole-1) 

Ba 

Tafel  slope  (V) 

c 

concentration  (molem-3) 

D 

diffusion  coefficient  (m2  s-1) 

w 

binary  diffusion  coefficient  (m2  s-1) 

E 

cell  voltage  (V) 

F 

Forchheimer  coefficient 

Faraday’s  constant,  9.6485309  x 

104  Cmole-1 

g 

acceleration  due  to  gravity  (m  s-2) 

G 

Gibbs  function  (Jkg-1) 

h 

convective  heat  transfer  coefficient 

(W  m-2  K"1) 

H 

enthalpy  (Jkg-1) 

I 

current  (A) 

j 

current  density  (Am-2) 

Jo 

exchange  current  density  (Am-2) 

mass  flux  of  species  a  due  to  molecular  diffu¬ 
sion  (kg  m  2  s-1) 

K 

forward  reaction  rate  constant  (ms-1  or 
molem-2  s-1 ) 

krk 

relative  permeability  of  phase  k 

K 

permeability  (m2) 

Ke 

equilibrium  constant  for  the  acid-base  reaction 
in  the  polymer  electrolyte 

2 

water  content  of  the  polymer  electrolyte 

^ k 

correction  applied  to  property  for  the  porous 
media  structure 

/\ 

M 

molecular  weight  (kg  mole-1) 

Mkn 

unit  normal  of  the  interface  between  phases  k 
and  n 

N 

molar  flux  (molem-2  s-1) 

P 

pressure  (Pa) 

p 

production  of  species  a  (molem-3  s-1) 

pH2(J 

Gat 

vapor  pressure  of  water  (Pa) 

q 

heat  flux  due  to  conduction  and  species  diffu¬ 
sion,  W  m-2 

Greact 

heat  of  reaction  (W  m-3) 

f react 

heat  of  reaction  (W  m-2) 

r 

radius,  m;  interaction  parameter  for  CO  ad¬ 
sorption  (J  mole-1) 

n 

reaction  rate  (molem-2  s-1) 

m 

universal  gas  constant,  8.31451  J  mole-1  K-1 

Sk 

saturation  of  phase  k 

S 

entropy  (Jkg-1  K-1) 

^energy 

source  term  in  the  conservation  of  energy  equa¬ 
tion  (W  m-3) 

^mom 

source  term  in  the  conservation  of  momentum 
equation  (N  m-3) 

^species 

source  term  in  the  conservation  of  species 
equation  (kgm-3  s-1) 

t 

time  (s) 

T 

temperature  (K) 

u 

velocity  (ms-1) 

Ukn 

velocity  of  the  interface  between  phases  k  and 
n  (ms-1) 

V 

volume  (m3) 

w 

relative  phase  velocity  (ms-1) 

W 

mass  flux  (kgm-2  s-1) 

xk 

mole  fraction  of  species  a 

Za 

charge  of  species  a 

Greek  letters 

p 

symmetry  factor  for  the  CO  adsorp¬ 
tion/desorption  reaction 

activity  coefficient  of  species  a 

,k 

interfacial  source  term  for  energy  conservation 
(J  m-3  s-1) 

r?,k 

interfacial  source  term  for  momentum  conser¬ 
vation  (N  m-3  s-1) 

interfacial  source  term  for  mass  conservation 
(kgm-3  s-1 ) 

r^a 

1  S  ,k 

interfacial  source  term  for  species  conservation 
(kgm-3  s-1) 

€k 

volume  fraction  of  phase  k 

T] 

overpotential  (V) 

ec 

contact  angle 

P 

fraction  of  reaction  sites  covered  by  species  a 

Ka,k 

electrical  conductivity  (S  m-1) 

k 

thermal  conductivity  (W  m-1  K-1) 

viscosity  (kgm-1  s-1) 

01 

P 

electro-chemical  energy  of  species  a 
(J  mole-1) 

P 

density  (kgm-3) 

a 

surface  tension  (N  m-1) 

rk 

viscous  stress  (N  m-2) 

0 

potential  (V) 

0 

void  fraction 

,  v a 

Mk 

mass  fraction  of  species  a 

Subscripts 

a 

anode 

c 

cathode;  capillary 

cell 

cell 

e 

polymer  electrolyte 

g 

gas 

i 

interface 

1 

liquid 

m 

mixture 

ref 

reference  value 
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rev  reversible 

s  solid 

Superscripts 
a  anode 

c  cathode 

e  polymer  electrolyte 

eff  effective  value 

Overbars 

per  unit  volume  (m-3) 
equilibrium  value 
correction  for  dispersion 
molar  units 

Operators 

( )  total- volume  average 

()*  phase- volume  average 


mathematical  models  in  the  published  literature,  with  varying 
degrees  of  simplifications. 

The  hydrogen  oxidation  and  oxygen  reduction  reactions  in 
the  catalyst  layers  were  modeled  by  Bernardi  and  Verbrugge 

[1] .  The  model  of  Bernardi  and  Verbrugge  [1]  assumed  that 
the  catalyst  layer  void  regions  were  filled  with  polymer  elec¬ 
trolyte;  other  models  allow  for  a  void  region  filled  with  gas 

[2]  or  a  combination  of  gas,  liquid  and  polymer  electrolyte 

[3] .  If  the  hydrogen  fuel  contains  CO,  the  hydrogen  oxidation 
reaction  is  inhibited,  severely  decreasing  the  efficiency  and 
power  output  of  the  PEM  fuel  cell.  This  CO-poisoning  of  the 
anode  was  modeled  by  Wang  and  Savinell  [4],  using  an  em¬ 
pirical  correlation  for  coverage  of  CO  on  the  anode  catalyst 
surface.  Springer  et  al.  [5]  addressed  CO-poisoning  by  us¬ 
ing  detailed  reaction  kinetics  for  the  adsorption,  desorption, 
and  electro-oxidation  of  hydrogen  and  carbon  monoxide.  The 
most  effective  method  of  mitigating  CO-poisoning  is  the  in¬ 
troduction  of  oxygen,  with  concentrations  ranging  from  1 
to  4%,  into  the  anode  fuel  stream;  this  is  referred  to  as  O2- 
bleeding.  The  adsorption,  desorption,  electro-oxidation,  and 
heterogeneous  oxidation  of  CO,  hydrogen,  and  oxygen  were 
incorporated  into  a  PEM  fuel  cell  model  by  Baschuk  and  Li 
[6]. 

Since  the  polymer  electrolyte  is  permeable  to  water,  pro¬ 
ton  migration  is  strongly  coupled  with  water  transport.  Ver¬ 
brugge  and  Hill  [7]  developed  a  mathematical  model  of  the 
proton  and  water  migration  in  the  pores  of  a  fully  humidified 
polymer  electrolyte  membrane.  However,  the  protonic  con¬ 
ductivity  is  a  function  of  the  membrane  hydration;  as  well, 
water  transport  can  be  driven  by  water  concentration  differ¬ 
entials.  Therefore,  Springer  et  al.  [8]  developed  a  PEM  fuel 
cell  model  in  which  the  protonic  conductivity  was  a  func¬ 
tion  of  membrane  hydration.  Water  transport  in  the  mem¬ 
brane  by  the  mechanisms  of  electro-osmotic  drag  and  dif¬ 


fusion  were  also  included  within  the  model.  However,  the 
model  of  Springer  et  al.  [8]  was  one-dimensional,  neglected 
pressure-driven  water  transport  and  relied  on  an  empirical 
correlation  for  the  conductivity  of  the  membrane.  Hence, 
Thampan  et  al.  [9]  used  the  Dusty  Gas  Model  and  acid- 
base  equilibrium  to  formulate  a  model  describing  proton 
conduction  and  water  transport  in  the  polymer  electrolyte 
membrane. 

Both  gas  and  liquid  phases  exist  within  the  PEM  fuel  cell, 
and  the  mass  transport  within  each  phase  affects  PEM  fuel 
cell  performance.  The  gaseous  reactants  in  the  PEM  fuel  cell 
must  travel  from  the  gas  flow  channels,  through  the  porous 
electrode  backing  layer,  and  into  the  catalyst  layers.  Gurau 
et  al.  [10]  developed  a  two-dimensional  model  that  incorpo¬ 
rated  the  reactant  flow  in  the  gas  flow  channels,  electrode 
backing  layers,  and  catalyst  layers  of  a  PEM  fuel  cell.  The 
catalyst  layer  formulation  was  similar  to  that  of  Bernardi  and 
Verbrugge  [1],  and  the  polymer  electrolyte  layer  formulation 
was  similar  to  that  of  Verbrugge  and  Hill  [7].  Um  et  al.  [11] 
developed  a  model  similar  to  that  of  Gurau  et  al.  [10],  but  al¬ 
lowed  for  the  simulation  of  the  transient  operation  of  a  PEM 
fuel  cell.  Siegel  et  al.  [12]  also  developed  a  two-dimensional 
model,  but  considered  the  catalyst  layer  void  regions  to  be 
composed  of  gas  and  polymer  electrolyte  membrane.  Addi¬ 
tionally,  proton  migration  and  water  transport  in  the  polymer 
electrolyte  membrane  layer  was  incorporated  through  a  for¬ 
mulation  similar  to  Springer  et  al.  [8].  Complex  gas  flow 
channel  designs  require  a  three-dimensional  model  in  order 
to  accurately  simulate  the  mass  transport.  Thus,  the  models  of 
Gurau  et  al.  [10]  and  Um  et  al.  [1 1]  were  extended  into  three 
dimensions  by  Zhou  and  Liu  [13]  and  Um  and  Wang  [14], 
respectively.  Shimpalee  et  al.  [15]  used  the  commercial  CFD 
software  FLUENT  to  model  a  PEM  fuel  cell;  the  catalyst 
layers  and  polymer  electrolyte  membrane  were  incorporated 
as  boundary  conditions.  Heat  transfer  can  be  significant  in  a 
PEM  fuel  cell,  especially  for  cells  operating  within  a  stack; 
thus  Berning  et  al.  [16]  incorporated  heat  transfer  into  a  three- 
dimensional  PEM  fuel  cell  model.  The  catalyst  layers  were 
considered  to  be  surfaces  and  modeled  with  boundary  condi¬ 
tions. 

Water  is  produced  in  the  cathode  catalyst  layer  as  a  product 
of  oxygen  reduction.  Due  to  the  polymer  electrolyte  mem¬ 
brane  hydration  requirements,  the  reactant  streams  typically 
enter  the  PEM  fuel  cell  fully  humidified;  hence  the  water  pro¬ 
duced  in  the  cathode  catalyst  layer  must  exit  the  PEM  fuel  cell 
in  liquid  form.  The  two-phase  flow  in  the  cathode  gas  flow 
channel  and  electrode  backing  layer  was  considered  by  Wang 
et  al.  [17].  Using  Darcy’s  law  and  the  Multi-phase  Mixture 
Model  (MMM)  [18],  liquid  water  transport  in  the  electrode 
backing  layer  was  expressed  as  a  function  of  the  capillary 
pressure  and  gravity.  He  et  al.  [19]  used  a  two-fluid  analysis  to 
model  the  liquid  water  transport  in  the  cathode  catalyst  layer; 
the  liquid  water  velocity  was  assumed  to  be  proportional  to 
the  gas  velocity  and  liquid  saturation  gradient.  Shimpalee  et 
al.  [20]  added  liquid  water  transport  to  the  model  of  Shim¬ 
palee  et  al.  [15]  by  treating  liquid  water  as  a  diffusing  species. 
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You  and  Liu  [21]  used  a  mixture  formulation,  similar  to  [18], 
to  add  liquid  water  transport  to  the  model  of  Gurau  et  al.  [  10] . 
Ferng  et  al.  [22]  modeled  the  two-phase  flow  in  the  cathode 
gas  flow  channel  and  electrode  backing  layer  by  assuming 
that  the  gas  and  liquid  phases  had  the  same  velocity.  The  wet¬ 
ting  nature  of  the  pores  has  a  significant  effect  on  the  capillary 
pressure,  and  thus  on  liquid  water  transport.  Hence,  Stockie 
[23]  simulated  liquid  water  transport  in  the  electrode  back¬ 
ing  layer  for  both  hydrophilic  and  hydrophobic  wettabilities. 
In  addition,  the  water  contact  angle  within  the  pores  and  the 
porosity  also  affect  the  capillary  pressure,  and,  correspond¬ 
ingly,  the  liquid  water  transport.  Therefore,  Pasaogullari  and 
Wang  [24,25]  modeled  liquid  water  transport  in  the  electrode 
backing  layer  for  different  contact  angles,  porosities,  and 
wettabilities. 

The  three  key  processes  in  a  PEM  fuel  cell,  the  electro¬ 
chemical  reactions,  proton  migration,  and  mass  transport, 
are  strongly  coupled.  However,  the  aforementioned  models, 
Refs.  [1-25],  either  utilize  a  formulation  that  does  not  model 
all  three  key  processes  in  sufficient  detail,  or  do  not  incor¬ 
porate  the  entire  PEM  fuel  cell  within  the  model.  Currently, 
including  all  three  key  processes  in  a  detailed  manner  is  diffi¬ 
cult,  since  no  general  formulation  exists  for  a  PEM  fuel  cell: 
governing  equations  are  currently  derived  for  a  specific  layer 
or  process  within  the  PEM  fuel  cell.  Therefore,  the  objective 
of  this  paper  is  to  present  a  comprehensive,  mathematical 
model  from  which  other  future  modeling  simulation  studies 
can  be  carried  out. 

The  physical  problem  of  the  PEM  fuel  cell  is  discussed, 
with  an  emphasis  placed  on  the  processes  occurring  within  the 
cell  and  the  structure  and  nature  of  the  flow  in  each  layer.  The 
input  and  output  parameters  for  the  model  are  chosen  to  cor¬ 
respond  with  the  dependent  and  independent  parameters  of 
a  PEM  fuel  cell  experimental  investigation.  Within  the  PEM 
fuel  cell,  several  phases  co-exist:  gas,  liquid,  and  solid.  The 
interactions  between  the  gas,  liquid  and  solid  phases  within 
the  PEM  fuel  cell  are  included  using  the  technique  of  volume 
averaging.  The  volume-averaged  conservation  equations  ap¬ 
ply  throughout  the  PEM  fuel  cell,  with  certain  terms  having 
unique  forms  in  each  layer  of  the  PEM  fuel  cell.  Since  the 
gas  and  liquid  phases  can  be  discontinuous  within  the  PEM 
fuel  cell,  a  mixture  approach  is  used  whereby  the  gas  and 
liquid  are  considered  as  a  single,  continuous  pseudo-fluid,  or 
mixture;  conservation  equations  are  then  developed  for  the 
pseudo-fluid. 

2.  Physical  problem 

A  PEM  fuel  cell  consists  of  several  components,  and  sev¬ 
eral  processes  occur  within  each  component,  as  illustrated  in 
Fig.  1.  The  reactant  gas  streams  enter  the  fuel  cell  through 
the  gas  flow  channels,  which  are  grooved  into  the  bipolar 
plates.  For  a  single  PEM  fuel  cell,  the  bipolar  plates  are  also 
referred  to  as  flow  distribution  plates.  The  bipolar  plates  are 
typically  made  of  graphite;  the  gas  flow  channels  are  ap- 


x 


h2,  h2o,  co2,  o2,  n2,  h2o 

co,n2,  o2 

Fig.  1 .  The  processes  occurring  within  a  PEM  fuel  cell.  The  PEM  fuel  cell  is 
composed  of  the  (a)  bipolar  plate,  (b)  gas  flow  channel,  (c)  electrode  backing 
layer,  (d)  catalyst  layer,  and  (e)  polymer  electrolyte  layer. 

proximately  1  mm  thick,  while  the  entire  bipolar  plate  has  a 
thickness  of  approximately  2  mm  [26].  The  main  reactant  on 
the  anode  side  of  the  PEM  fuel  cell  is  hydrogen.  However, 
if  the  hydrogen  is  formed  from  a  hydrocarbon  source,  car¬ 
bon  dioxide  and  carbon  monoxide  can  also  be  present.  The 
presence  of  carbon  monoxide,  in  concentrations  of  greater 
than  2  ppm,  inhibits  the  oxidation  of  hydrogen  in  the  anode 
catalyst  layer  (CO-poisoning)  [27].  Introducing  oxygen,  in 
concentrations  of  1-4%,  into  the  anode  gas  stream  can  mit¬ 
igate  CO-poisoning  (02 -bleeding  [28]);  hence  oxygen  and 
nitrogen  can  be  present  in  the  anode  gas  flow  channels.  Oxy¬ 
gen  is  the  reactant  on  the  cathode  side  of  the  cell,  and  nitro¬ 
gen  can  be  present  if  air  is  used  as  the  oxidant.  The  polymer 
electrolyte  requires  humidification;  hence  both  the  anode  and 
cathode  streams  are  generally  fully  humidified. 

In  order  to  reach  the  catalyst  layers,  the  reactants  must  pass 
through  the  porous,  electrode  backing  layers.  Constructed  of 
carbon  paper  or  cloth,  the  electrode  backing  layers  are  ap¬ 
proximately  200  jmm  in  thickness  [26] .  The  electrode  backing 
layers  also  allow  liquid  water  to  exit  the  catalyst  layers  and 
enter  the  gas  flow  channels. 

The  conversion  of  the  chemical  energy  of  the  reactants  into 
electrical  energy,  heat,  and  liquid  water  occurs  in  the  catalyst 
layers,  which  have  a  thickness  of  approximately  5  jam  [26]. 
The  catalyst  layers  are  also  porous;  reactant  gas,  liquid  water, 
and  polymer  electrolyte  occupy  the  void  space  and  the  solid 
matrix  consists  of  carbon-supported  platinum  catalyst.  If  the 
fuel  is  CO-free,  the  overall  reaction  (Ra)  in  the  anode  catalyst 
layer  is  hydrogen  oxidation: 

H2 — >  2H+  +  2e~.  (1) 

The  polymer  electrolyte  is  electronically  insulative;  hence  the 
electron  produced  by  the  above  hydrogen  oxidation  must  pass 
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through  the  electrode  backing  layer,  bipolar  plate,  and  the 
external  load.  The  proton  in  reaction  (1)  is  transferred  through 
the  polymer  electrolyte  and  participates  in  the  reduction  of 
oxygen  in  the  cathode  catalyst  layer: 

2H+  +  2e“  +  lo2  — >  H20(i).  (2) 

The  combination  of  reactions  (1)  and  (2)  is  the  overall  PEM 
fuel  cell  reaction: 


H2  +  -O2 


>  H20(i)  +  heat  +  electrical  energy. 


The  heat  produced  by  the  electro-chemical  reactions  must 
exit  the  PEM  fuel  cell.  The  heat  can  be  transported  through 
conduction  in  the  solid  phase  of  the  electrode  backing  layers 
and  bipolar  plates,  or  through  convection  and  conduction  in 
the  fluid  phases. 

The  polymer  electrolyte  layer  consists  of  a  sulfonated  fluo- 
ropolymer,  which  is  similar  to  Teflon.  The  fluorinated  carbon 
chains  terminate  in  SO3H  groups,  which  when  in  contact  with 
water,  dissociate  into  SO3-  and  H30+  ions.  The  presence  of 
the  hydronium  ions  allow  the  sulfonated  fluoropolymer  to  act 
as  an  electrolyte;  the  protons  produced  in  the  anode  catalyst 
layer  can  be  conducted  through  the  membrane  in  order  to 
react  in  the  cathode  catalyst  layer.  The  electrolyte  is  perme¬ 
able  to  reactant  gases  and  water,  although  the  permeability 
to  reactant  gases  is  low  and  reactant  cross-over  is  small.  The 
protonic  conductivity  of  the  membrane  decreases  if  the  water 
content  of  the  membrane  is  low;  hence,  external  humidifica¬ 
tion  from  the  reactant  streams  is  necessary  in  order  to  prevent 
electrolyte  dehydration.  The  thickness  of  the  polymer  elec¬ 
trolyte  layer  ranges  from  50  to  250  [xm  [29]  and  the  most 
popular  polymer  electrolyte  employed  for  PEM  fuel  cells  is 
Nation. 

Each  layer  of  the  PEM  fuel  cell  has  a  unique  structure,  as 
shown  in  Fig.  2.  The  anode  and  cathode  reactants  typically 
enter  the  gas  flow  channels  fully  humidified.  As  a  result, 
the  water  produced  by  the  PEM  fuel  cell  must  be  entrained 


within  the  gas  flow  channels  as  small  liquid  droplets. 
Two-phase,  porous  media  flow  occurs  within  the  electrode 
backing  layers,  with  the  liquid  water  being  driven  by  capil¬ 
lary  pressure.  The  electrode  backing  layer  is  manufactured 
such  that  there  is  a  mix  of  hydrophobic  and  hydrophilic 
pores. 

The  transport  of  water,  protons,  and  reactant  gas  in  the 
polymer  electrolyte  can  be  explained  by  the  ionic  cluster 
model  of  Gierke  and  Hsu  [30].  In  the  presence  of  water,  the 
Nation  membrane  is  composed  of  three  regions  [31].  Liquid 
water  and  hydronium  ions  are  contained  within  clusters  that 
have  a  diameter  of  approximately  4  nm  [30].  Channels  inter¬ 
connect  the  clusters,  with  the  channels  having  a  diameter  of 
approximately  1  nm  [30].  Both  the  channels  and  the  clusters 
are  lined  with  immobile  sulfonate  ions.  Between  the  clus¬ 
ters  and  the  rigid,  hydrophobic  backbone  of  the  polymer,  ex¬ 
ists  the  amorphous  part  of  the  perfluorinated  backbone.  This 
hydrophobic  region  is  permeable  to  gases,  thus  transport  of 
reactant  gas  occurs  in  this  region. 

The  flow  characteristics  in  the  catalyst  layer  are  a  combi¬ 
nation  of  the  flow  in  the  electrode  backing  layers  and  poly¬ 
mer  electrolyte  layer.  Two-phase,  porous  media  flow  occurs 
in  the  electrolyte-free  void  space,  while  liquid  water  and  re¬ 
actant  gas  are  also  transported  within  the  electrolyte.  Electro¬ 
chemical  reactions  occur  at  the  catalyst/electrolyte  interface, 
thus  the  reactant  gases  are  consumed  and,  in  the  cathode  cat¬ 
alyst  layer,  liquid  water  is  produced. 


3.  Model  input  and  output 

The  input  and  output  parameters  of  the  PEM  fuel  cell 
model  are  chosen  with  regard  to  an  experimental  investigation 
into  PEM  fuel  cell  performance.  The  performance  of  a  PEM 
fuel  cell  can  be  characterized  by  a  voltage  versus  current  plot. 
When  testing  a  PEM  fuel  cell,  several  operating  conditions 
can  be  altered,  such  as 


Electrode  Backing 
Layer 


Polymer  Electrolyte 
Layer 


Fig.  2.  The  nature  of  the  flow  in  the  layers  of  a  PEM  fuel  cell. 
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•  the  inlet  concentration,  temperature,  and  flow  rate  of  the 

reactants; 

•  the  cell  temperature; 

•  the  electrical  load;  and 

•  either  the  inlet  or  outlet  pressure. 

The  temperature  variation  within  the  fuel  cell  is  small. 
Therefore,  during  operation,  the  temperature  is  measured  at 
one  location  within  the  cell  and  referred  to  as  the  cell  tempera¬ 
ture.  This  cell  temperature  is  controlled  by  heating  or  cooling 
the  external  surfaces  of  the  PEM  fuel  cell.  The  electrical  load 
on  the  PEM  fuel  cell  is  specified  by  setting  the  voltage  of  the 
PEM  fuel  cell  or  the  current  drawn  from  the  cell.  The  inlet  and 
outlet  pressures  are  related  through  the  pressure  drop  in  the 
gas  flow  channels;  hence,  specifying  both  pressures  as  input 
conditions  is  redundant.  If  the  PEM  fuel  cell  is  not  pressur¬ 
ized,  then  the  gas  flow  channels  exhaust  to  the  atmosphere; 
thus  the  outlet  pressure  is  specified.  During  pressurized  op¬ 
eration,  the  inlet  pressure  is  generally  specified  and  regulated 
with  a  back-pressure  control  valve  that  is  located  at  the  outlet 
of  the  gas  flow  channels. 

Therefore,  one  set  of  input  parameters  are  the  operating 
conditions,  and  these  become  the  external  boundary  condi¬ 
tions  for  the  PEM  fuel  cell  model.  The  location  of  the  external 
boundary,  for  a  two-dimensional  geometry,  is  shown  as  the 
dark,  solid  line  in  Fig.  3.  The  PEM  fuel  cell  in  Fig.  3  has  a 
thickness  of  Tl  and  a  height  of  Xl.  The  interfaces  between 
the  various  layers  of  the  PEM  fuel  cell  are  also  shown  in 
the  figure.  The  bipolar  plate/gas  flow  channel,  gas  flow  chan¬ 
nel/electrode  backing  layer,  electrode  backing  layer/catalyst 
layer,  and  catalyst  layer/polymer  electrolyte  layer  interfaces 
in  the  anode  side  of  the  PEM  fuel  cell  are  denoted  by  p/fc, 

^fc/eb’  ^eb/cl  ’  an<^  ^cl/e’  respectively.  Similar  notation  is  used 
to  denote  the  interfaces  on  the  cathode  side,  except  that  the 
superscript  “a”  is  replaced  by  “c”. 

The  external  boundary  conditions  can  be  classified  into 
three  groups:  the  x  surfaces  excluding  the  gas  flow  channels 
(i.e.,  0  <  y  <  Fbap/fc,  Fa/eb  <  y  <  Yfc/cb,  and  Fbcp/fc  <  y  < 
Tl);  the  y  surfaces;  and  the  x  surfaces  of  the  gas  flow  chan¬ 


Fig.  3.  The  boundaries  where  the  boundary  conditions  for  the  PEM  fuel  cell 
model  are  imposed.  The  dark  solid  line  denotes  the  external  boundary,  while 
the  dashed  line  is  the  location  of  the  internal  boundary. 


nels  (i.e.,  Tap/fc  <  y  <  Tfa/eb  and  F£/eb  <  y  <  Tbcp/fc).  The 
v  surfaces,  excluding  the  gas  flow  channels,  are  considered  as 
solid,  insulated  walls.  Therefore,  the  velocity,  heat  transfer, 
and  species  flux  are  zero: 


0  <  v  <  Ta 
u  -  y  —  2bp/fc 

yfc/eb  <  y  <  ffc/eb 

ybP/fc  <  y  <  yL 
x  =  0,  XL 


where  uk  is  the  velocity  of  phase  k,  the  diffusive  flux  of 
species  a  in  phase  k,  Js  the  current  density  in  the  solid  phase, 
and  qk  the  heat  flux  due  to  conduction  and  species  diffusion. 
The  unit  vector  in  the  x-direction  is  expressed  as  i.  The  bound¬ 
ary  condition  on  the  heat  flux  is  an  approximation,  since  the 
v  surfaces  of  a  single  PEM  fuel  cell  are  generally  exposed 
to  the  ambient  environment.  However,  the  aspect  ratio  of  the 
cell,  Tl/Xl,  is  generally  much  less  than  1.  Additionally,  the 
temperature  of  a  single  PEM  fuel  cell  is  controlled  by  heat¬ 
ing  or  cooling  the  y  surface  of  the  bipolar  plate.  Therefore, 
the  heat  lost  or  gained  from  the  x  surfaces  can  be  considered 
negligible,  compared  to  the  heat  transfer  at  the  y  surfaces  of 
the  bipolar  plates. 

Since  the  y  surfaces  are  heated  or  cooled,  boundary  con¬ 
ditions  on  the  heat  flux  must  be  specified.  The  heat  transfer 
at  the  bipolar  plates  can  be  specified  as  either  a  fixed  sur¬ 
face  temperature,  a  specified  heat  flux,  or  a  convective  heat 
transfer  condition: 


Ts  =  specified 
qs  =  specified 
qs-f±h(Ts  -  T^) 


=  0 


for 


0  <  x  <  XL 
y  =  0,  Tl 


where  is  the  ambient  air  temperature,  j  the  unit  vector  in 
the  y  direction,  and  h  the  convective  heat  transfer  coefficient. 
The  cell  voltage  can  be  expressed  as 

E  —  ^rev  Ocelli  (4) 


where  EYQV  is  the  reversible  cell  potential  and  rjce.\\  the 
cell  voltage  loss,  or  overpotential.  The  cell  overpotential 
includes  the  activation  overpotential  associated  with  the 
electro-chemical  reactions,  the  concentration  overpotential 
caused  by  mass  transport,  and  the  ohmic  losses  incurred  by 
electron  and  proton  transport  within  the  PEM  fuel  cell.  The 
reversible  cell  potential  can  be  found  by  assuming  that  the 
overall  fuel  cell  reaction  is  in  equilibrium: 


✓V  /V 

where  AG  is  the  change  in  Gibbs  energy  and  AS  the  change 
in  entropy  for  the  overall  fuel  cell  reaction.  The  inlet  tem¬ 
perature  and  pressure  of  the  reactants  are  denoted  by  T  and 
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P,  respectively;  the  universal  gas  constant  and  the  Faraday 
constant  are  given  as  Wl  and  respectively.  The  subscript 
“ref”  denotes  a  reference  value.  Using  the  standard  values  for 

/V  /V 

AG  and  AS,  the  reversible  cell  potential  can  be  calculated  as 
Er  =  1.229  -  0.85  x  10“3(T  -  298.15) 

+  4.31  X  10_5r  ln[(PH2)(Po2)1/2L 

where  T  is  in  K  and  P  is  in  atmospheres. 

The  absolute  value  of  the  potential  in  the  solid  phase  is  ar¬ 
bitrary,  and  depends  on  the  choice  of  reference  electrode  [32]. 
In  this  model,  the  reference  electrode  is  chosen  to  consist  of 
the  solid  phase  at  the  open  circuit  conditions.  Therefore,  the 
difference  between  the  potentials  on  the  anode  and  cathode 
bipolar  plates  becomes  the  cell  overpotential  [33].  Assum¬ 
ing  that  the  potential  on  the  surface  of  the  bipolar  plates  is 
constant,  the  boundary  conditions  for  the  cell  potential  are 

^s|y=0  =  7?cell>  ^sl  y=Y]^  =  0’  (6) 

where  cathode  potential,  <Ps\y=YL,  is  arbitrarily  defined  as 
zero.  Hence,  in  the  PEM  fuel  cell  mathematical  model,  the 
cell  voltage  is  calculated  by  specifying  the  overpotential.  In 
this  approach,  the  electric  double  layers  in  the  electrolyte  near 
the  two  catalyst  layers  have  been  neglected.  The  double  lay¬ 
ers  provide  the  driving  potential  needed  for  proton  transport 
through  the  electrolyte,  but  do  not  contribute  to  the  overall 
overpotential. 

The  anode  and  cathode  gas  flow  channels  each  have  one 
inlet  and  one  outlet.  The  inlets  can  be  either  on  the  same  side, 
which  is  referred  to  as  a  co-flow  arrangement,  or  on  opposite 
sides,  which  is  referred  to  as  a  counter-flow  arrangement. 
At  the  inlet,  the  velocity,  species  mass  fractions,  tempera¬ 
tures,  and  volume  fractions  are  specified  for  the  gas  and  liquid 
phases: 


Uk 

oc 

the  inlet  1 

cor 

rf 

h 

►  =  specified 

for 

r?o  a-  tfl 

& 

A  A 

Vi  Vi 

A  A 

v-  v- 

o1  a* 

where  €k  is  the  volume  fraction  and  the  mass  fraction  of 
species  a  within  phase  k  is  denoted  by  co%. 

As  mentioned  previously,  the  pressure  can  be  specified  at 
either  the  inlet  or  the  outlet;  here,  it  is  assumed  that  the  cell 
exhausts  to  the  atmosphere  and  the  pressure  is  specified  at 
the  outlet.  The  profiles  of  velocity,  mass  fraction,  and  tem¬ 
perature  are  not  known  at  the  outlet;  however,  if  the  gas  flow 
channels  are  extended  beyond  the  PEM  fuel  cell,  the  flow  be¬ 
comes  fully  developed.  Thus,  the  outlet  boundary  conditions 
can  be  expressed  as 


P  =  specified 


d_ 

dx 


X— Xl+^oo 


for 


I  the  outlet 

%/fc  <  <  Ifc/eb 

*?c/eb  <  y  <  Ibp/fc 


9 


where  is  the  extension  of  the  gas  flow  channels  that  re¬ 
sults  in  fully  developed  flow. 

The  external  boundary  conditions  coincide  with  the  oper¬ 
ating  conditions  of  a  PEM  fuel  cell.  However,  internal  bound¬ 
ary  conditions,  shown  as  the  dashed  line  in  Fig.  3,  are  also 
required,  due  to  phase  discontinuities.  The  H30+  ions  are 
only  present  in  the  polymer  electrolyte  and  catalyst  layers; 
hence  a  no-flux  boundary  condition  must  be  applied  at  the 
electrode  backing/catalyst  layer  interfaces: 

wf3°+  -j  =  0  for  v  =  Feab/Cl  and  Feb/cl, 


where  W^3°  is  the  mass  flux  of  hydronium  ions.  Likewise, 
electrons  cannot  be  transported  in  the  polymer  electrolyte 
layer;  hence  a  no-flux  boundary  condition  is  also  necessary 
at  the  catalyst  layer/polymer  electrolyte  membrane  layer  in¬ 
terfaces: 

Js  J  =  o  for  y=  K‘l|/e  and  Fcc1/e. 

Since  the  bipolar  plate  is  impermeable,  the  fluid  phases  dis¬ 
appear  at  the  bipolar  plate/gas  flow  channel  interface.  There¬ 
fore,  the  appropriate  boundary  conditions  are  the  no-slip  and 
no- temperature  jump  conditions: 

for  y  -^bp/fc  Tbp/fc. 

The  performance  can  also  be  influenced  by  the  manufac¬ 
ture  of  the  PEM  fuel  cell;  this  effect  is  included  within  the 
model  by  input  parameters  that  are  referred  to  as  design  pa¬ 
rameters.  These  include  such  parameters  as  the  overall  di¬ 
mensions  of  the  PEM  fuel  cell  and  the  thickness  of  each  layer. 
The  choice  of  materials  for  the  bipolar  plates  influences  the 
heat  transfer  and  electrical  conductivity  properties  of  the  cell. 
For  the  porous  electrode  backing  layer,  the  porosity  and  per¬ 
meability  are  important  input  parameters,  as  well  as  the  nature 
of  the  pores:  the  degree  of  hydrophobic  and  hydrophilic  be¬ 
havior  influences  the  capillary  pressure.  The  fraction  of  the 
void  space  occupied  by  polymer  electrolyte  in  the  catalyst 
layer  is  also  an  important  design  parameter.  In  addition  to  the 
design  parameters,  the  model  requires  values  for  the  phys¬ 
ical  properties  of  the  reactant  gas  and  liquid  water,  such  as 
density,  viscosity,  diffusivity,  and  thermal  conductivity. 

The  mathematical  model  analyzes  the  transport  of  mass, 
momentum,  species,  and  energy  within  the  gas,  liquid,  and 
solid  phases  of  the  PEM  fuel  cell.  Therefore,  the  velocity, 
pressure,  concentration,  potential,  and  temperature  profiles 
within  the  PEM  fuel  cell  are  the  model  output  parameters. 
From  these  profiles,  the  mass  flux  of  each  species  and  the 
current  density  can  be  determined.  This  information  is  sig¬ 
nificant  because,  due  to  the  small  geometry  and  operating 
environment,  direct  in  situ  measurement  of  these  quantities 
within  a  PEM  fuel  cell  is  difficult,  if  not  impossible.  Since 
the  cell  voltage  is  an  input  parameter,  the  current  output  of 
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the  PEM  fuel  cell  is  determined  with 

I=(  JsjdA,  (7) 

^  A>p 

where  Abp  is  the  surface  area  of  the  bipolar  plate  and  Js  the 
current  density,  which  changes  over  the  active  cell  surface. 

Therefore,  the  mathematical  model  has  two  inputs:  the 
operating  conditions  and  the  design  parameters.  The  operat¬ 
ing  conditions  correspond  to  the  independent  variables  during 
the  experimental  investigation  of  a  PEM  fuel  cell.  The  design 
parameters  depend  on  the  manufacture  of  the  PEM  fuel  cell. 
Profiles  of  velocity,  pressure,  concentration,  potential,  and 
temperature  are  determined  by  the  model  as  output  parame¬ 
ters.  Since  the  cell  voltage  is  specified  as  an  input  parameter, 
the  performance  of  the  PEM  fuel  cell  is  quantified  by  the 
corresponding  cell  current,  or  current  density  distribution. 


4.  Formulation 


species  diffusion.  The  enthalpy  of  the  phase  does  not  include 
any  contribution  from  charged  species;  hence 

Hk  =  £  o/kH«,  (12) 

0't^Q'± 

where  H ®  is  the  enthalpy  of  species  a  in  phase  k  and  a± 
represents  a  charged  species.  If  the  polymer  electrolyte  is 
free  of  contaminants,  then  the  only  charged  species  are  the 
H30+  ions  in  the  liquid  phase  of  the  polymer  electrolyte.  The 
energy  of  a  charged  species  is  potential-dependent;  therefore, 
it  is  represented  by  the  E term  in  the  conservation  of  energy 
equation: 

Ef  =  ~wf  -V Hf .  (13) 


Using  the  Gibbs  equation,  the  enthalpy  of  a  charged  species 
can  be  expressed  in  terms  of  the  electro-chemical  energy 

4-  /v  dz 

(ji®  )  and  the  partial  molar  entropy  (S%  )  [32]: 


(14) 


Three  phases  co-exist  within  the  PEM  fuel  cell:  gas,  liquid, 
and  solid.  The  conservation  of  mass,  momentum,  species,  and 
energy  within  the  phases  are  expressed  with  [34] 

^  +  V  •  (pkuk)  =  0,  (8) 


The  electro-chemical  energy  can  be  expressed  as 
V£f  =  m,TkV  In  af  + 

where  za t  is  the  charge  and  af  is  the  activity  [35]: 


(15) 

(16) 


— C okuk )  +  V  •  ( pkukuk )  +  VPk  -  V  •  rk  -  pkg  -  b  =  0, 
ot 

(9) 

^(Pkcoak)  +  V  •  (pkcoakuk  +  J?)  =  0,  (10) 

ot 

A PkHk )  +  V  •  (pkHkuk)  =  -V-qk  +  Ef,  (11) 

ot 

where  Eqs.  (8)— ( 11)  are  the  conservation  equations  for  mass, 
momentum,  species,  and  energy,  respectively.  The  subscript 
k  denotes  the  phase,  with  k  =  g,  1  and  s  representing  the  gas, 
liquid,  and  solid  phases,  respectively.  Although  Eqs.  (8)— (1 1) 
apply  to  all  three  phases,  the  immobility  of  the  solid  phase 
results  in  only  the  conservation  of  species  and  energy  being 
relevant  to  the  solid  phase. 

The  density  and  velocity  of  the  phase  are  denoted  by  pk  and 
Uk ,  respectively.  In  the  conservation  of  momentum,  Eq.  (9), 
momentum  change  due  to  pressure  (Pk),  viscous  stress  (r^), 
gravity  (g),  and  a  body  force  (b)  are  included.  It  is  assumed 
that  no  reactions  take  place  within  either  the  gas  or  liquid 
phases;  thus,  the  conservation  of  species  a  within  phase  k 
is  represented  by  Eq.  (10).  The  mass  fraction  of  species  a 
within  phase  k  is  denoted  by  go while  represents  the 
mass  flux  of  species  a  due  to  molecular  diffusion.  Energy 
transfer  due  to  pressure  work  and  viscous  dissipation  are  both 
neglected  in  the  conservation  of  energy,  which  is  expressed 
as  Eq.  (11).  The  enthalpy  of  phase  k  is  denoted  by  Hk,  and 
qk  represents  the  transport  of  energy  due  to  conduction  and 


The  activity  coefficient  is  y®  and  x®  is  the  mole  fraction.  The 
definition  of  activity  given  in  Eq.  (16)  is  valid  for  both  charged 
and  uncharged  species  and  for  an  ideal  gas  or  a  dilute  solution, 
the  activity  coefficient  is  unity.  Therefore,  the  energy  of  the 
charged  species  can  be  expressed  as 

±  Wa±  ±  ± 

Eak  = - 4-  ■  mrkV  In aak  +  za±^V<Pk  +  V(TkSak  )]. 

Ma 

(17) 

For  the  solid  phase,  the  charged  species  are  electrons  and 
only  the  potential  gradient  term  is  relevant;  thus  Eq.  (17) 
represents  Joule  heating  in  the  solid  phase. 

The  viscous  stress,  diffusional  mass  flux,  and  the  heat  flux 
due  to  conduction  and  species  diffusion  are  expressed  as  [34] 


rk  =  fik[^uk  +  (Va*)^, 

(18) 

— c“V  In  a\  -  mTkV<Pk 

=  T  ,  kk  (Jkjka  jf, 

MaMpWa-P'k  k  k 

(19) 

qk  =  -XkVTk  +  JfH«Jka. 

(20) 

a 


In  the  expression  for  viscous  stress,  f  represents  the  transpose 
and  ilk  is  the  viscosity  of  phase  k;  it  is  assumed  that  each  phase 
is  a  Newtonian  fluid  and  that  the  dilation  effect  on  the  viscous 
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stress  tensor  is  negligible.  The  Generalized  Stefan-Maxwell 
equations  are  used  to  describe  the  diffusive  mass  flux,  with  the 
temperature  and  pressure  effects  neglected.  Since  electrons 
are  the  only  mobile  species  in  the  solid  phase,  the  general¬ 
ized  Stefan-Maxwell  equation  reduces  to  Ohm’s  law  in  the 
solid  phase.  The  molar  concentration  is  ck  and  the  molecular 

/K  /K 

weights  of  the  phase  and  species  are  denoted  by  M \  and  Ma, 
respectively.  In  Eq.  (20),  Xk  is  the  thermal  conductivity. 

The  transfer  of  mass,  momentum,  species,  and  energy 
from  adjacent  phases  is  absent  from  the  conservation  equa¬ 
tions  presented  in  Eqs.  (8)— (1 1).  In  order  to  account  for  these 
processes,  the  conservation  equations  must  be  integrated  over 
a  representative  volume: 


(conservation  equation) 


1 


(conservation  equation)  dV, 


(21) 


where  VY  is  the  representative  volume.  Using  the  transport 
and  averaging  theorems  [36],  the  volume-averaged  conserva¬ 
tion  equations  can  be  expressed  in  terms  of  volume-averaged 
quantities,  such  as  density  and  velocity,  and  interfacial  source 
terms.  For  quantities  such  as  density  and  velocity,  the  phase- 
volume,  or  intrinsic,  average  is  defined  as 


<**>*  =  IT  (  ^dy’ 
yk  Jvk 


(22) 


where  Vk  is  the  volume  of  phase  k  within  the  representative 
volume  and  44  is  a  quantity  within  phase  k ,  such  as  density 
or  velocity.  The  phase- volume  average  of  Eq.  (22)  and  the 
total- volume  average  of  Eq.  (21)  are  related  with 


(Vk)  =  Vk 

(%}*  Vr 


(23) 


where  <4  is  the  volume  fraction  of  phase  k  within  the  repre¬ 
sentative  volume. 

The  general,  volume-averaging  procedure  for  two-phase, 
gas/liquid  flows  is  described  by  Ishii  [37],  while  Slattery  [38] 
presents  the  volume-averaging  procedure  for  porous  media 
flows.  Additionally,  volume- averaging  is  applied  to  the  mod¬ 
eling  of  batteries  and  the  porous  region  of  a  PEM  fuel  cell  by 
Wang  et  al.  [39].  After  the  volume-averaging  procedure,  the 
conservation  equations  become 


Utk(Pk)*)  +  V  •  (€k{pk)*{uk)*)  =  rM,k , 

ot 


(24) 


U^kiPk)* {uk)*)  +  V  •  (€k(pk)*(uk)*(uk)*) 

ot 

+  Vfe(ft)*)  -  <P*)*V(e*)  -  V  •  (€k(Tk)*) 


%-(€k(Pk)*(Hk)*)  +  V  •  (€k(pk)*(Hk)*(uk)*) 
ot 

=  -V  •  (€k{qk)*)  +  €k(Eff  +  rEik,  (27) 

where  r  represents  the  interfacial  source  terms.  The  volume- 

I  * 

averaged  form  of  the  energy  of  the  charged  species,  (Ek  )  , 
is  assumed  to  be  the  same  as  Eq.  (17),  except  that  the  bulk 
values  are  replaced  by  volume-averaged  values.  The  volume- 
averaged  forms  of  the  viscous  stress  tensor  ((4)*)  and  heat 
flux  ((q)k)  are 

(xk)*  =  i4lf[V(uk)*  +  (V<«jfc>*)+],  (28) 

(q)t  =  -*f  V<7*>*  +  (HZ)*(Jk)*,  (29) 


where  the  effective  viscosity  and  thermal  conductivity  are 
denoted  by  /xjrff  and  A^ff,  respectively.  The  effective  values 
of  viscosity  and  thermal  conductivity  depend  on  the  nature 
of  the  multiphase  flow;  however,  they  have  the  general  form 
of  [40] 

=  (1  +  L%)Vk  +  (30) 

where  Lf  is  the  correction  for  the  porous  media  structure  and 

V-/ 

44  the  correction  for  dispersion  and  hydrodynamic  effects. 

Because  of  the  small  pore  diameters,  especially  within  the 
polymer  electrolyte,  molecular  interactions  between  species 
in  the  fluid  phases  and  the  solid  phase  can  be  significant. 
These  interactions  can  be  incorporated  within  the  diffusional 
mass  flux  term  by  considering  the  solid  phase  as  a  diffus¬ 
ing  species  with  a  velocity  of  zero  [41].  Thus,  the  volume- 
averaged  form  of  the  generalized  Stefan-Maxwell  equations 
become 


<c?>*V  In  (4)*- 

Mk 


Za(cak)*d W 


(Tk) 
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V(0e> 
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=  E 
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MaMsDeit  R, 

ffa  a  P  a-p,k 
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M  Z)eff 

1V1°(Ua-s,k 


(31) 


where  is  the  total  mass  flux  of  species  a  and  D^_  k  repre¬ 
sents  the  interactions  between  the  solid  phase  and  the  species 
within  the  fluid  phases.  In  regions  of  the  PEM  fuel  cell  where 
the  solid  phase  does  not  exist,  or  does  not  have  an  apprecia¬ 
ble  effect  on  species  mass  flux,  DQ^_s  k  — >  00  and  Eq.  (31) 
reduces  to  the  standard  Stefan-Maxwell  equation:  Eq.  (19). 
An  equation  for  the  diffusive  mass  flux  can  be  obtained  by 
inverting  Eq.  (31): 


-€k(Pk)*g-€k(b)*  =  rFtk, 

jt(€k{pk)*(afk)*)  +  V  •  (ek(pk)*(a>%)*(uk)* 

+  ek(jkr)  =  rlk. 


(25) 
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(32) 


where  D^fk  is  the  overall  diffusion  coefficient  of  species  a  in 

phase  k  and  K^k  is  the  electrical  conductivity  for  species  a. 
These  terms,  as  well  as  Tk,  are  defined  as 


1  1 

Deff  Deff  ^  Deff  ’ 

^ a—f>,k  ^  a— s,k 

eff  _  (Za)2(ct)*^Dfk 

<a'k  m(Tk)*o  -  (xak)*y 


(33) 


(34) 


Eqs.  (24)  and  (27)  represent  the  conservation  of  mass, 
momentum,  species,  and  energy  in  the  gas,  liquid  and  solid 
phases  of  the  PEM  fuel  cell.  In  order  to  solve  these  equations, 
constitutive  equations  are  required  for  several  parameters. 
Some  constitutive  equations  apply  throughout  the  PEM  fuel 
cell;  these  are  equations  relating  the  temperatures  of  the  fluid 
phases,  the  concentration  of  water  in  the  gas  phase,  and  the 
interfacial  source  term  between  the  fluid  and  solid  phases 
in  the  conservation  of  momentum  equation.  Local  thermal 
equilibrium  is  assumed  between  the  gas  and  liquid  phases; 
therefore,  the  temperatures  are  equal: 

(Ei)*  =  (Eg)*.  (41) 

Also,  if  the  liquid  phase  is  present,  the  partial  pressure  of 
water  within  the  gas  phase  is  equal  to  the  vapor  pressure  in 
the  water;  therefore, 


r f  =  l  + 


91n(^>* 

9  In  (*£)*' 


(35) 


Mathematically,  the  interfacial  source  terms  in  Eqs.  (24)- 
(27)  are  expressed  as 


PkiMk  Mkn)  ’  M kn  dA, 


PkUk(uk  Ukn)  '  Mkn  dA 


Pk^kn  dA 


(36) 


Tk  •  njcn  dA, 


Pk^k(Mk  Ukn)  ’  Mkn  dA, 


PkHk(uk  Ukn)  '  4lkn  dA 


qk  •  nkn  dA, 


(37) 

(38) 


(39) 


where  represents  the  summation  over  all  adjacent  n 
phases,  Akn  the  interfacial  area  of  phases  k  and  n,  and  ukn  the 
velocity  of  the  interface  of  k  and  n  within  the  representative 
volume.  The  pressure  deviation,  Pk,  is  given  by  [36] 

h  =  Pk-(Pk )*•  (40) 


PH20 

( JC?2°)*  =  -2**- 

g  (Pg>* 


(42) 


The  momentum  interfacial  source  term  between  the  fluid  and 
solid  phases  can  be  modeled  with  the  Darcy-Forchheimer 
terms  [42]: 


Em,£  =  - 


(€k)2r-k 

Kkrk 


(»*>*  - 


\J  Kkrk 


\(uk)*\(uk)*,  (43) 


* 


where  K  is  the  absolute  permeability  of  the  porous  media, 
krk  the  relative  permeability,  and  F  the  Forchheimer  term.  If 
no  solid  phase  exists,  then  K  =  oo  and  F  =  0.  The  absolute 
permeability  and  the  Forchheimer  term  depend  on  the  struc¬ 
ture  of  the  porous  media,  while  the  relative  permeabilities  are 
functions  of  ek.  The  F  term  is  not  expected  to  be  significant 
in  the  polymer  electrolyte,  due  to  the  low  fluid  velocities. 
However,  it  will  be  significant  in  the  electrode  backing  layer 
in  the  region  adjacent  to  the  gas  flow  channels  [42]. 

The  form  for  most  of  the  constitutive  equations  required 
for  closure  of  Eqs.  (24)— (27)  are  location-specific,  with  dif¬ 
ferent  forms  being  required  in  the  polymer  electrolyte  layer, 
electrode  backing  layer,  gas  flow  channel,  and  catalyst  layer. 
Specifically,  expressions  for  the 


•  Volume  fraction; 

•  Body  forces; 

•  Pressure  difference  between  the  phases; 

•  Reaction  kinetics  and  corresponding  interfacial  source 
term  for  the  conservation  of  species  equation;  and 

•  Interfacial  source  term  for  the  energy  equation 


Physically,  Fy[,k  represents  the  mass  entering  phase  k  from 
all  adjacent  phases.  The  momentum  transfer  from  the  adja¬ 
cent  phases  to  phase  k  due  to  mass  transfer,  pressure,  and 
viscous  forces  is  represented  by  Fptk.  Species  transfer  from 
other  phases  is  expressed  as  Fgk,  while  is  the  interfacial 
energy  transfer  due  to  mass  transfer  and  molecular  conduc¬ 
tion. 


are  required.  However,  not  all  of  the  above-mentioned  param¬ 
eters  are  required  in  each  layer.  The  body  force  is  only  sig¬ 
nificant  in  the  polymer  electrolyte  and  catalyst  layers,  while 
only  a  relationship  between  the  gas  and  liquid  pressures  is 
required  within  the  gas  flow  channels.  The  relevant  constitu¬ 
tive  equations  for  each  layer  are  examined  in  the  following 
sections. 
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4.1.  Formulation  for  the  polymer  electrolyte 


The  gas,  water,  and  proton  transport  properties  of  the  poly¬ 
mer  electrolyte  are  functions  of  the  membrane  water  content 
(2),  which  is  defined  as 


moles  of  liquid  water 
moles  of  SO3H 


(44) 


As  the  electrolyte  absorbs  more  water,  the  value  of  2  in¬ 
creases.  The  transport  of  water  and  ions  in  the  electrolyte 
can  only  occur  if  the  water  content  is  above  a  minimum  level 
(2min)  [5].  Hence,  the  percolation  model  should  be  applied 
to  the  effective  diffusion  coefficients  within  the  electrolyte 
[9],  whereby  of_jU  =  0  if  S<  /mm- 

The  volume  fraction  of  the  membrane  available  for  water 
and  proton  transport  becomes  [9] 


2 

(t>s/t>H2o)  +  ^ 


(45) 


(a)  Hydrophilic  pore  (b)  Hydrophobic  pore 


Fig.  4.  Cylindrical  pore  having  (a)  hydrophilic  and  (b)  hydrophobic  charac¬ 
teristics. 


Using  a  mass  balance,  the  mole  fractions  of  hydronium  and 
liquid  water  in  the  pores  of  the  electrolyte  become  a  function 
of  the  membrane  hydration  [9] : 


,  h3o+x*  1 

<  V  )  = 


(1  +  2)  -  7(1  +  2e)2  -  42(1  -  K~l) 


2(1  -  K”1) 


(49) 


/V 

where  Vs  is  the  partial  molar  volume  of  the  polymer  and 

✓V 

Vh7o  the  partial  molar  volume  of  liquid  water.  The  perme¬ 
ability  of  the  electrolyte  to  both  gas  and  liquid  increases  as 
the  water  content  increases  [31];  this  implies  that  the  relative 
permeability  and  the  volume  fraction  of  the  gas  phase  are 
proportional  to  the  membrane  water  content: 


where  Ke  is  the  equilibrium  constant  for  reaction  (48). 

The  small  pore  size  makes  the  assumption  of  local  ther¬ 
modynamic  equilibrium  between  the  solid  and  fluid  phases 
within  the  polymer  electrolyte  layer  reasonable.  Therefore, 
within  the  polymer  electrolyte  layer, 

(7g)*  =  (71)*  =  <rs>*.  (50) 


krk  —  ft2, 


4.2.  Formulation  for  the  electrode  backing  layer 


where  the  constant  of  proportionality  is  expressed  as  f.  The 
percolation  model  should  also  be  applied  to  the  relative  per¬ 
meability  of  the  liquid  water;  hence,  kri  =  0  if  2  <  2min. 

An  electrical  potential  gradient  exists  within  the  polymer 
electrolyte,  and  this  gradient  exerts  a  body  force  on  the  wa- 
ter/hydronium  ion  mixture.  Therefore,  within  the  liquid  phase 
of  the  polymer  electrolyte,  the  body  force  is  defined  as  [7] 

(*i>*  =  -<cfl3°+>: Vv<0e>*,  (46) 

where  c[*3  °  is  the  molar  concentration  of  H3  0+  in  the  liquid 
phase  and  <£e  the  potential  in  the  polymer  electrolyte.  The 
body  force  of  Eq.  (46)  is  responsible  for  the  electro-osmotic 
drag  effect  in  the  polymer  electrolyte. 

The  relationship  between  the  pressures  in  the  gas  and  liq¬ 
uid  phases  are  unknown;  due  to  the  small  geometry  of  the 
pores,  direct  measurement  is  impossible.  However,  since  the 
area  in  which  the  gas  transport  occurs  is  flexible,  it  will  be 
assumed  that  the  pressures  in  the  gas  and  liquid  phases  are 
equal.  Hence 

(rPg)*  =  (Pi)*.  (47) 

The  concentrations  of  H30+  and  liquid  water  are  related 
through  the  acid-base  equilibrium  reaction  occurring  within 
the  electrolyte  pores  [9]: 

SO3H  +  H2O  ^  H30+  +  SO3".  (48) 


In  the  electrode  backing  layer,  the  volume  fraction  of  the 
solid  phase  is  a  design  parameter.  Therefore,  it  is  convenient 
to  define  the  saturation  of  the  gas  or  liquid  phase  as 

Sk  =  (51) 

<P 

where  the  void  fraction  (0)  is  the  volume  available  to  the  fluid 
phases,  with  0  =  1  —  6S. 

The  pressures  in  the  gas  and  liquid  phases  are  related  with 
the  capillary  pressure: 


(52) 


Both  the  relative  permeability  (krk)  and  the  capillary  pressure 
are  a  function  of  the  liquid  saturation  in  the  porous  media.  In 
a  cylindrical  pore,  the  difference  between  the  gas  and  liquid 
pressure  is  a  function  of  the  surface  tension,  contact  angle, 
and  pore  radius  [43]: 


2(7  cos  0C 

p  — _ _ 

rc  — 

r 


(53) 


where  0  is  the  surface  tension  between  the  gas  and  liquid 
phases,  0C  is  the  contact  angle,  and  r  is  the  pore  radius.  If  the 
pore  material  is  hydrophilic,  the  contact  angle  is  less  than  90° 
and  Pc  is  positive;  hydrophobic  materials  have  a  negative  Pc 
and  the  water  contact  angle  is  greater  than  90°,  as  illustrated 
in  Fig.  4. 
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A  porous  media  can  be  considered  to  be  composed  of  a 
large  distribution  of  pores  with  different  dimensions,  with 
the  water  content  of  each  pore  depending  on  the  applied 
capillary  pressure.  The  variation  of  capillary  pressure  with 
water  content  in  porous  media  is  illustrated  in  Fig.  5.  For 
hydrophilic  porous  media,  the  capillary  pressure  is  always 
positive,  as  shown  in  Fig.  5(a),  while  Fig.  5(b)  illustrates 
that  hydrophobic  porous  media  are  characterized  by  negative 
capillary  pressures.  Mixed- wettability  porous  media,  shown 
in  Fig.  5(c),  can  exhibit  both  positive  and  negative  capillary 
pressures  [44]. 

All  three  capillary  pressure  curves  have  two  similar  fea¬ 
tures.  The  first  is  the  existence  of  hysteresis,  with  different 
curves  being  produced  if  water  is  entering  (wetting)  or  leav¬ 
ing  (drying)  the  porous  media;  thus  the  capillary  pressure 
depends  on  both  liquid  saturation  and  history.  The  hysteresis 
is  due  to  the  difference  in  contact  angle  for  advancing  and 
receding  water  [45]. 

The  second  common  feature  between  the  two  capillary 
pressure  curves  is  that  the  capillary  pressure  asymptotically 
approaches  infinity  at  a  finite,  minimum  value  of  saturation, 
and  approaches  negative  infinity  for  a  maximum  value  of  wa¬ 
ter  saturation.  This  can  be  interpreted  as  there  being  a  residual 
saturation  of  water,  %,  and  a  residual  saturation  of  gas,  sgr, 
in  the  porous  media.  The  residual  saturation  represents  the 
saturation  at  which  the  fluid  loses  the  capability  to  move  as  a 
bulk  phase  in  response  to  a  hydraulic  gradient,  resulting  in  the 
transport  of  the  residual  fluid  being  dominated  by  advective- 
diffusive  transport  as  a  dispersed  phase  in  the  fluid  with  the 
larger  volume  fraction  [46] . 


Attempts  have  been  made  to  correlate  capillary  pressure 
curves,  with  a  comprehensive  representation  being  [46] 


Pnw 


/\y  =  22,(7  COS  (0cw) 


5 


(54) 


where  Pnw  and  Pw  are  the  pressures  of  the  non- wetting  and 
wetting  fluids,  respectively,  Z  is  a  correction  factor  to  account 
for  the  change  in  contact  angle  due  to  roughness,  0C w  is  the 
intrinsic  contact  angle  of  the  wetting  fluid  to  the  solid,  and 
C(s^)  is  a  function  of  the  effective  saturation  of  the  wetting 
fluid.  The  effective  saturation  is  defined  as 


Vw  -  s Wr 
1  —  %Wr  —  AVr  ’ 


(55) 


where  ^NWr  and  swr  are  the  residual  saturations  of  the  non¬ 
wetting  and  wetting  fluids,  respectively.  Therefore,  if  one 
capillary  pressure  curve  is  available,  it  can  be  applied  to  dif¬ 
ferent  porous  media,  using  corrections  for  the  porosity,  per¬ 
meability  and  roughness.  However,  the  function  C  is  not  fully 
universal,  and  depends  on  several  other  parameters,  such  as 
pore  structure. 

Predictive  models  for  the  relative  permeability  were 
developed  from  conceptual  models  of  flow  in  capillary  tubes 
combined  with  models  of  pore-size  distribution  [47].  The 
common  conceptual  models  are  the  Burdine  and  Mualem 
functions.  In  addition  to  the  conceptual  models,  several 
porous  media  modeling  studies  have  successfully  used  a  sim¬ 
ple  power  law  relationship  for  the  relative  permeability  [40] : 

krg  =  ( sff ,  kr 1  =  <f)n.  (56) 


Fig.  5.  Variation  of  capillary  pressure  with  liquid  saturation  for  (a)  hydrophilic,  (b)  hydrophobic,  and  (c)  mixed- wettability  porous  media  [44].  The  pressure 
head  difference  is  equal  to  Pc/(p\g),  where  Pc  is  the  capillary  pressure  and  g  is  the  acceleration  due  to  gravity. 
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The  temperature  of  the  solid  and  fluid  phases  may  be  un¬ 
equal;  thus,  an  expression  for  the  interfacial  heat  transfer  be¬ 
tween  the  solid  and  fluid  phases  is  required.  One  method  is 
to  use  a  convective  heat  transfer  coefficient  [48]: 

rE,k  =  hksAks((Ts )*  -  (7*)*),  (57) 

s 

where  hks  is  the  convection  coefficient  and  A^s  is  the  inter¬ 
facial  area  per  unit  volume.  Correlations  for  the  convection 
coefficient  are  available  from  the  published  literature  [40]. 

4.3.  Formulation  for  the  gas  flow  channel 

No  solid  phase  exists  in  the  gas  flow  channels;  hence, 
no  interfacial  source  terms  are  required.  Since  the  water  is 
assumed  to  be  in  suspended  droplet  form,  the  relationship 
between  the  pressures  can  be  expressed  as 

(Pg)*  -  (P\)*  =  ~  — ,  (58) 

r 

where  o  is  the  surface  tension  and  r  is  the  characteristic  radius 
of  the  droplet. 

4.4.  Formulation  for  the  catalyst  layer 

The  catalyst  layer  is  a  porous  media,  similar  to  the  elec¬ 
trode  backing  layer;  however,  the  structure  is  more  compli¬ 
cated  due  to  the  fact  that  polymer  electrolyte  exists  within  the 
pore  region.  Fig.  6  shows  that  the  gas  and  liquid  phases  are 
contained  within  two  structures:  the  void  region  of  the  cata¬ 
lyst  layer  and  the  polymer  electrolyte.  Therefore,  five  sets  of 
conservation  equations  are  applied  in  the  catalyst  layer,  with 


1 .  k  =  g  referring  to  the  gas  phase  in  the  void  region; 

2.  k  =  1  referring  to  the  liquid  phase  in  the  void  region; 

3.  k  =  g,  e  referring  to  the  gas  phase  within  the  polymer 
electrolyte; 

4.  k  =  1,  e  referring  to  the  liquid  phase  within  the  polymer 
electrolyte; 

5.  k  =  s  referring  to  the  carbon  support  and  catalyst. 

The  conservation  equations  for  k  =  g,  1,  and  s  are  similar 
to  the  corresponding  equations  in  the  electrode  backing  layer, 
while  the  k  =  g,  e  and  1,  e  cases  are  similar  to  the  equations 
in  the  polymer  electrolyte  layer.  The  volume  fractions  can  be 
expressed  as 

€g  =  (j)Sg,  6i  =  (j>S\ ,  6g,e  =  #e£g,  fl,e  =  <pSte\, 

where  0  =  1  —  es  and  e®  and  are  the  volume  fractions  of 
the  electrolyte  that  contain  gas  and  liquid  phase,  respectively. 
The  pressure  relationship  between  the  gas  and  liquid  phases  in 
the  void  region  are  the  same  as  in  the  electrode  backing  layer, 
while  the  pressure  in  the  polymer  electrolyte  is  considered  to 
be  the  result  of  the  gas  and  liquid  phase  pressures: 

(PgW  =  (P\,W  =  ~  ( Pg ) *  +  (59) 

The  water  content  of  the  polymer  electrolyte  is  determined 
by  the  water  content  of  the  gas  and  liquid  phases  within  the 
void  region  of  the  catalyst  layer: 

2=  (60) 
1  Sq 

where  Sg  and  2j  are  the  hydrations  when  the  membrane  is  in 
contact  with  a  gas  or  liquid  phase,  respectively.  The  value  of 
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Fig.  6.  The  structure  of  the  catalyst  layer. 
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2j  can  be  considered  as  a  function  of  temperature  only  [49], 
but  the  value  of  Sg  is  a  function  of  both  temperature  and  the 
concentration  of  water  vapor  in  the  gas  phase.  The  value  of 
can  be  found  using  a  Flory-Huggins  model  [50]. 
Electro-chemical  and  heterogeneous  reactions  occur  on 
the  interface  between  the  polymer  electrolyte  and  the  solid 
phases  in  the  catalyst  layers.  The  interfacial  source  terms  in 
the  conservation  of  species  equation  represent  the  rate  of  pro¬ 
duction  of  each  species;  this  rate  of  production  is  determined 
by  the  reaction  kinetics.  In  the  anode  catalyst  layer,  the  fol¬ 
lowing  reactions  are  considered  [6] : 

1.  Hydrogen  adsorption,  desorption  and  electro-oxidation; 

2.  Carbon  monoxide  adsorption,  desorption  and  electro¬ 
oxidation; 

3.  Heterogeneous  oxidation  of  carbon  monoxide  by  oxygen; 

4.  Heterogeneous  oxidation  of  hydrogen  by  oxygen. 


Hence,  the  reactions  can  be  expressed  as 
H2  +  2M  ^  2(H-M),  (61) 

2(H-M)  +  2H20  ^  2M  +  2H30+  +  2e“,  (62) 

CO  +  M  ^  (CO-M),  (63) 

(CO-M)  +  3H20  ^  M  +  C02  +  2H30+  +  2e“,  (64) 

02  +  2M  ^  2(0-M),  (65) 

(CO-M)  +  (O-M)  ->  C02  +  2M,  (66) 

(O-M)  +  2(H-M)  -*  H20  +  3M,  (67) 


where  reactions  (61)  and  (62)  represent  the  adsorption  and 
electro-oxidation  of  hydrogen,  assumed  to  follow  Langmuir 
and  Butler- Volmer  kinetics,  respectively;  reactions  (63) 
and  (64)  are  the  adsorption  and  electro-oxidation  of  carbon 
monoxide,  assumed  to  follow  Temkin  and  Butler- Volmer 
kinetics,  respectively;  and  reactions  (65)-(67)  denote  the 
heterogeneous  oxidation  of  hydrogen  and  carbon  monoxide, 
assumed  to  follow  Langmuir-Hinshelwood  kinetics.  The 
reaction  rates,  in  units  of  mole  m-2  s-1,  are 


7?H  _  kU  rrH2r6)Mx2  _  jlH  r6)Hx2x 

7Vl,ads  —  )  ^adsWs  )  1’ 
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{^sh)z-cxr(oji, 
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a,  ox 


"OX 


-  *„*  uc!(?'ri 


ox 
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(71) 

(72) 

(73) 

(74) 


where  k f  denotes  the  forward  rate  reaction  constant,  bf  de¬ 
notes  the  ratio  of  backward  to  forward  reaction  rate  constant, 
and  is  the  exchange  current  density.  The  fraction  of  the 
platinum  reaction  sites  in  the  solid  phase  covered  by  species 
a  is  denoted  by  with^  =  1  —  6^  —  0jr°  —  6®  represent¬ 
ing  the  fraction  of  free  reaction  sites.  The  interface  tempera¬ 
ture  is  represented  by  7)  and  the  concentrations  are  assumed 
to  be  equivalent  to  the  volume-averaged  concentrations.  An 
overbar  over  the  concentration  or  coverage,  represents 
the  value  at  equilibrium,  or  when  the  reaction  rates  are  zero. 

The  overpotential  of  the  anode  catalyst  layer  is  denoted  as 
?7a  and  is  defined  as  [32] 


?7a  —  <Z>s  f7a, 


(75) 


where  the  value  of  U'  is  taken  as  zero.  In  the  carbon  monox- 

a 

ide  adsorption/desorption  reaction  rate,  Eq.  (70),  p  is  a  sym¬ 
metry  factor  that  has  a  value  between  zero  and  one,  and 
r  is  an  interaction  parameter  that  represents  the  effect  of 
lateral-interaction  on  the  adsorption/desorption  process  [51]. 
In  Eqs.  (69)  and  (71),  B a  is  the  Tafel  slope  for  species  a ; 
the  Tafel  slopes  for  the  forward  and  backward  directions  are 
considered  equivalent,  with 


B"  =  Bf 


(76) 


The  production  of  reaction  intermediates  (0^,  0jr°,  and 
0p)  are  assumed  to  be  in  steady  state;  thus, 


a  —  27^aads 
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CO 

s  =  Vco  -  V CO  -  Vco~°  =  0 
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(79) 


Using  Eqs.  (77)-(79),  the  production  of  the  reactants  and 
products  can  be  expressed  in  terms  of  ft? ox,  ft? ?x,  ft?0? 

and  nfL0- 
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,H 


--1Z 
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-n 
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The  major  reaction  occurring  in  the  cathode  catalyst  layer 
is  oxygen  reduction: 

02  +  4H30+  +  4e“  =  6H20.  (80) 


The  rate  of  reaction  is  governed  by  Volmer/Erdey-Gruz  ki¬ 
netics  [52]: 
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(81) 


where  is  the  exchange  current  density  for  oxygen  reduc¬ 
tion.  The  Tafel  slope  is  denoted  by  B®2  and  is  given  by 


mr 

dd 


(82) 


Eq.  (82)  results  in  a  Tafel  slope  of  68  mV/decade  at  a  tem¬ 
perature  of  70  °C,  which  agrees  with  the  experimental  value 
of  70  mV/decade  from  Parthasarathy  et  al.  [53]. 

Gases  can  dissolve  into  the  polymer  electrolyte  and  be 
transported  between  the  anode  and  cathode;  however  trans¬ 
port  of  gases  in  the  polymer  electrolyte  is  slow.  Since  CO  is 
present  in  only  a  ppm  amount,  it  is  unlikely  that  a  significant 
amount  could  be  transported  across  the  polymer  electrolyte 
membrane  layer  and  react  in  the  cathode  catalyst  layer.  How¬ 
ever,  cross-over  of  hydrogen  could  be  possible  and  as  a  re¬ 
sult,  the  heterogeneous  oxidation  of  hydrogen  by  oxygen  is 
included  in  the  cathode  catalyst  layer.  The  chemical  reaction 
and  reaction  rate  are  represented  in  the  same  manner  as  in 
the  anode  catalyst  layer;  hence,  the  net  production  of  each 
species  in  the  cathode  catalyst  layer  is 
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Using  the  reaction  kinetics,  the  interfacial  source  term  in 
the  conservation  of  species  equation  is  related  to  the  iP 
terms  through 

rs“*  =  AnMaPa,  (83) 


(84) 


S 

where  An  is  the  reactive  surface  area  per  volume  in  the  cata- 

Q  1  /V 

lyst  layer.  The  units  of  rgk  are  kg  m  s  ,  while  r®k  is  in 

molar  units  of  mole  m~3  s-1. 

The  reactions  occurring  within  the  catalyst  layers  are 
exothermic;  thus,  the  heat  of  reaction,  Greact,  is  added  to  the 
conservation  of  energy  in  the  solid  phase  [48]: 

/ 

Greact  —  AntheacU  (85) 


<7re act  =  -  ^  H“] 
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—  greact,  ir  T  greact, rev  T  greact,  het, 


(86) 


where  greact, ir  and  greact. rev  are  the  irreversible  and  reversible 
heat  produced  by  the  electro-chemical  reactions.  The  heat 
produced  by  the  heterogeneous  chemical  reactions  is  ex¬ 
pressed  as  greact,  het-  The  irreversible  heat  production  is  a  func¬ 
tion  of  the  overpotential: 


greact,  ir  —  ‘ 
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The  reversible  heat  production  is  a  function  of  the  entropy 
change  of  the  electro-chemical  reactions: 
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The  change  in  entropy  of  the  reactions,  A  Sr,  are 
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Finally,  the  heat  produced  by  the  heterogeneous  oxidation  of 
hydrogen  and  carbon  monoxide  is  expressed  as 
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where  the  enthalpy  changes  are 

A  £/H-0  _  fjH20  _  £yH2  _  1  lr02 

A7/r  —  Hl,e  Hg,e  2  6e’ 
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5.  Final  set  of  governing  equations 

The  previous  sections  outlined  the  governing  equations 
for  the  conservation  of  mass,  momentum,  species  and  energy 
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in  the  gas  and  liquid  phases,  as  well  as  the  conservation  of 
species  and  energy  in  the  solid  phase.  Because  the  gas  and  liq¬ 
uid  phases  can  be  non-continuous,  the  fluid  phases  are  com¬ 
bined  to  form  a  pseudo-fluid,  or  mixture.  Thus,  the  processes 
occurring  within  the  PEM  fuel  cell  can  be  described  by  two 
sets  of  equations,  with  one  set  applying  to  the  pseudo-fluid 
and  the  other  set  applying  to  the  solid  phase.  The  pseudo-fluid 
equations  apply  from  T£p/fc  <  y  <  y£p/fc,  while  the  equa¬ 
tions  for  the  solid  phase  apply  between  0  <  y  <  T^e  and 

^l/e  <  y  <  YL. 

The  properties  of  the  pseudo-fluid  are  defined  as 

Pm  =  Sg(Pg)*  +  Sl(pi)*  +  Se€eg{pg,e)*  +  5e«f<Pl,e)*>  (95) 

PmMm  =  Sg<Pg»g>*  +  Sl<Pl»l>*  +  ^e^(Pg,e)*(«g,e>* 

+  Se4(Pl,e»l,e>*,  (96) 

Pm  =  Sg(Pg)*  +  S^)*  +  Se6'(Pg,e>*  +  5e*f  (97) 

Pm«“  =  Sg(Pg)*  (C0ag)*  +  S,<p,)*«>*  +  See>g,e>’V“e>* 

+  ^4<Pl,e)*<<e>*-  (98) 

PmHm  =  Sg(pg)*(Hg)*  +  Si(pi)*  W*  +  Se6eg(pg)*  (Hg)* 

+  se^(pl)*(Hl)*,  (99) 

where  pm ,  Hm,  Pm  >  co ^  and  Hm  are  the  density,  velocity,  pres¬ 
sure,  mass  fraction  and  enthalpy  of  the  mixture,  respectively. 

Using  the  definitions  of  the  mixture  properties,  the  con¬ 
servation  of  mass,  momentum,  species,  and  energy  in  the 
pseudo-fluid  can  be  expressed  as 

—  (0Pm)  +  V  •  (0pmwm)  =  0,  (100) 
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is  equivalent  to  the  potential  in  the  liquid  phase  of  the  poly¬ 
mer  electrolyte  ((0e)*).  The  mixture  viscous  stress,  diffusion 
coefficient,  and  thermal  conductivity  are  given  with 
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where  the  summation  is  over  the  gas  and  liquid  phases. 

The  source  term  in  the  momentum  equation,  <Smom,m,  has 
three  components: 
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where  rm  is  the  viscous  stress,  D^fm  is  the  effective  diffusiv- 

ity  of  species  a ,  and  is  the  thermal  conductivity  in  the 
pseudo-fluid.  In  Eq.  (102),  the  potential  in  the  mixture  ( @m ) 


where  S^om  m  represents  momentum  transfer  due  to  inter¬ 
facial  pressure  differences,  5^^  m  represents  momentum 
transfer  due  to  the  solid  surface  and  the  electrical  body  force, 
and  S^om  m  represents  momentum  transfer  due  to  relative 
phase  motion.  The  permeability  and  Forchheimer  term  for 
the  mixture  are 
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The  relative  phase  velocities,  w,  are  defined  as 


(113) 
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The  source  term  for  the  conservation  of  species  equation, 
iSSpecies,m>  can  also  be  decomposed  into  three  components: 
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(119) 


where  5^Lc-es  m  represents  the  production  or  consumption  of 

species  a  due  to  the  electro-chemical  reactions,  are 

the  non-Fickian  terms  of  the  Stefan-Maxwell  equations,  and 
<srel  .  m  is  due  to  relative  phase  motion.  While  the  source 

in 

terms  due  to  the  non-Fickian  terms  and  relative  phase  motion 


are  continuous,  the  production  and  consumption  of  reactants 
occurs  only  in  the  catalyst  layers.  Hence,  Species  m  on^ 
non-zero  in  the  catalyst  layers. 

The  source  term  in  the  conservation  of  energy  equation, 
^energy, m,  can  be  considered  as  the  sum  of  four  terms: 
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^energy.m  =  hmAs((Ts)*  -  Tm),  (121) 
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where  Energy  m  the  convective  heat  transfer  between  the 

solid  and  the  mixture,  <S^3e?gy,m  represents  heat  generation 

due  to  H30+  migration,  ^energy  m  due  t0  relative  phase 
motion,  and  <S^nergy  m  is  the  energy  stored  in  the  polymer 
backbone  of  the  polymer  electrolyte.  The  volume  fraction, 
density,  and  specific  heat  of  the  polymer  are  €eM,  pM,e,  and 
Cp >g,  respectively.  The  surface  area  of  the  solid  phase  is  de- 

s 

noted  by  As. 

The  source  terms  in  the  mixture  conservation  equations 
are  influenced  by  the  relative  phase  velocities;  expressions 
for  w  are  required  for  closure.  In  order  to  develop  expres¬ 
sions  for  the  relative  velocities,  assumptions  must  be  made 
on  the  nature  of  the  flow  in  each  component  of  the  PEM  fuel 
cell.  The  flow  in  the  gas  flow  channels  is  considered  to  be 
homogeneous;  thus,  wg  =  w\  =  0.  In  the  electrode  backing 
layers,  the  flow  of  liquid  water  is  expected  to  be  small.  Hence, 
the  relative  velocity  of  the  liquid  phase  can  be  determined  by 
assuming  that  motion  is  governed  by  Darcy’s  law: 

®1  =  -;^v  «fl>*  -  <P1> d25) 

(pSlJLl 

The  relative  velocity  of  the  gas  phase  can  be  determined  with 
the  relationship  of  Eq.  (115). 

In  a  similar  manner,  the  relative  velocity  of  the  liquid  phase 
in  the  polymer  electrolyte  layer  can  be  determined  by  assum¬ 
ing  that  the  flow  can  be  described  by  the  Schlogl  equation: 

=  _KJc rl£v((jp  }*  _  b  _  {  y  )  _  ^  (126) 

<q>l,e 

The  relative  velocity  of  the  gas  phase  within  the  polymer 
electrolyte  can  be  calculated  using  Eq.  (115). 

The  catalyst  layer  requires  values  for  w\,  ufi,e  and  wg?e, 
since  wg  can  be  calculated  with  Eq.  (115).  The  relative  ve¬ 
locity  of  the  liquid  in  the  catalyst  layer  void  space  can  be 
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found  with  Eq.  (125),  while  w\^e  is  given  by  Eq.  (126).  As¬ 
suming  that  the  flow  velocity  of  the  gas  phase  within  the 
polymer  electrolyte  is  small,  the  relative  velocity  of  the  gas 
phase  within  the  polymer  electrolyte  is 

Wg,e  =  _^MV((/>g,e>*  -  (ps,e)*g)  -  um.  (127) 

€g^  g,e 

For  the  solid  phase,  only  the  conservation  of  species  and 
energy  are  applicable.  The  flux  of  electrons  is  expressed  in 
terms  of  current  density,  using  Faraday’s  Law,  and  because 
electrons  are  the  only  mobile  species  in  the  solid  phase, 
the  Generalized  Stefan-Maxwell  equations  reduce  to  Ohm’s 
Law.  Hence,  the  conservation  of  species  and  energy  can  be 
simplified  to 

0  =  -V  •  (es<ffV(<Z>s)*)  +  <Sspecies,j!  (128) 

0  =  V  •  (eskf- V<rs)*)  +  Energy, „  (129) 

where  <SSpecies,s  and  <Senergy,s  are  source  terms.  The  source 
term  in  the  conservation  of  species  equation  includes  the  con¬ 
sumption  or  production  of  electrons,  while  the  source  term 
for  the  conservation  of  energy  equation  includes  the  Joule 
heating  term,  convective  heat  transfer  from  the  mixture,  and 
the  heat  of  the  reaction: 


S, 


species,  s 


(130) 


^energy,  s  —  € 


J, 


+  hmAs(Tm  -  <rs)*)  +  2 react-  (131) 


S 


The  source  term  for  the  conservation  of  species  equation  is 
only  non-zero  in  the  catalyst  layers;  the  heat  of  reaction  is 
also  non-zero  in  the  catalyst  layers. 

Therefore,  the  processes  occurring  within  a  PEM  fuel  cell 
can  be  modeled  with  two  sets  of  equations.  The  first  set, 
Eqs.  (100)-(103)  describe  the  transport  of  mass,  momentum, 
species,  and  energy  within  a  mixture  consisting  of  the  gas  and 
liquid  phases.  The  second  set,  Eqs.  (128)  and  (129),  represent 
the  migration  of  electrons  and  energy  transport  in  the  solid 
phase  of  the  PEM  fuel  cell.  If  there  are  N  uncharged  species 
within  the  gas  and  liquid  phases  of  the  PEM  fuel  cell,  the 
two  sets  of  equations  will  sum  to  a  total  of  N  +  5  equations. 
These  equations  can  be  solved  for  the 


•  mixture  velocity  («m)  and  pressure  (Pm), 

•  N  —  1  mixture  mass  fractions  of  the  uncharged  species 

«), 

•  potential  in  the  solid  phase  ((<£)*)  and  electrolyte  ( @m ), 
and 

•  mixture  enthalpy  (Hm)  and  solid  phase  temperature 

«rs}*). 


However,  the  source  terms  contained  within  the  equation 
sets  are  functions  of  quantities  within  the  individual  phases, 
or  phase  quantities.  As  well,  the  mixture  density  is  a  func¬ 
tion  of  the  liquid  saturation.  Thus,  the  solution  of  the  mixture 
conservation  equations  is  iterative,  and  the  mixture  quanti¬ 
ties  must  be  converted  into  phase  quantities  at  each  itera¬ 


tion.  The  first  step  in  the  conversion  between  the  mixture 
and  phase  quantities  is  the  determination  of  the  liquid  phase 
saturation  (s\).  This  is  accomplished  by  setting  a  =  H2O 
in  the  definition  of  the  mixture  mass  fraction,  Eq.  (98).  If 
liquid  water  is  present,  the  mass  fraction  of  the  water  va¬ 
por  can  be  determined  with  Eq.  (42)  and,  if  the  electrolyte 
phase  is  also  present,  the  mass  fraction  of  liquid  water  within 
the  electrolyte  can  be  calculated  by  the  acid-base  equilib¬ 
rium  reaction.  Since  sg  +  s\  =  1,  Eq.  (98)  can  be  solved 
for  s\. 

Once  the  liquid  saturation  is  determined,  the  mixture  den¬ 
sity  can  be  calculated  with  Eq.  (95).  Also,  the  pressure  of 
the  gas  and  liquid  phases  can  be  determined  with  the  mix¬ 
ture  pressure  and  the  capillary  pressure  functions.  With  the 
individual  phase  pressures  and  Eqs.  (125)-(127),  the  rela¬ 
tive  velocities  of  the  gas  and  liquid  phases  can  be  calculated, 
along  with  the  actual  velocities  of  the  gas  and  liquid  phases. 
The  mass  fractions  of  each  species  within  the  gas  and  liquid 
phases  can  be  determined  with  the  mixture  mass  fractions 
and  Eq.  (98).  The  temperatures  of  the  gas  and  liquid  phases 
are  assumed  equal;  hence,  using  equations  of  state  and  the 
definition  of  the  mixture  enthalpy,  Eq.  (99),  the  temperature 
of  the  gas  and  liquid  phases  can  be  found.  Using  the  liquid 
saturation,  the  relative  velocities  of  each  phase,  and  the  in¬ 
dividual  phase  values  of  pressure,  temperature,  and  species 
mass  fractions,  the  source  terms  in  the  mixture  and  solid  phase 
conservation  equations  can  be  evaluated.  Therefore,  new  val¬ 
ues  of  the  mixture  velocity,  mixture  pressure,  mixture  mass 
fractions,  potentials,  mixture  enthalpy,  and  solid  phase  tem¬ 
perature  can  be  evaluated. 

After  convergence  of  the  iterative  solution  process,  the 
mixture  quantities  can  be  converted  into  phase  quantities, 
allowing  for  the  display  and  analysis  of  the  velocity,  pres¬ 
sure,  mass  fractions,  and  temperatures  within  the  individual 
phases.  The  cell  voltage  is  specified  as  in  input  condition; 
hence,  the  cell  performance  can  be  quantified  by  the  cell  cur¬ 
rent,  which  is  calculated  using  Eq.  (7). 


6.  Conclusions 

A  mathematical,  PEM  fuel  cell  model  was  developed  that 
described  the  conservation  of  mass,  momentum,  species  and 
energy  within  the  bipolar  plates,  gas  flow  channels,  electrode 
backing,  catalyst,  and  polymer  electrolyte  layers.  The  model 
consisted  of  two  equation  sets.  The  first  equation  set  described 
the  transport  of  mass,  momentum,  species,  and  energy  within 
a  mixture  consisting  of  a  gas  and  liquid  phase.  Electron  mi¬ 
gration  and  energy  transport  in  the  solid  phase  of  the  PEM 
fuel  cell  were  considered  with  the  second  set  of  equations. 
The  governing  equations  for  the  mixture  and  solid  were  de¬ 
rived  using  volume-averaging  techniques. 

The  conservation  equations  for  the  mixture  included 
source  terms  that  represented  fluid/solid  interactions  and  ef¬ 
fects  due  to  relative  phase  motions.  Thus,  the  source  terms 
described  the  relevant  phenomena  that  was  significant  in  each 
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layer  of  the  PEM  fuel  cell.  Darcy-Forchheimer  drag  terms 
were  included  in  order  to  account  for  the  effect  of  the  dis¬ 
persed  solid  phase  on  the  momentum  of  the  fluid  phases 
within  the  electrode  backing,  catalyst,  and  polymer  elec¬ 
trolyte  layers.  Electro-osmotic  drag  in  the  electrolyte  was 
incorporated  with  an  electrical  body  force.  Also,  the  con¬ 
sumption  and  production  of  species  within  the  catalyst  layers 
were  included  as  a  source  term  in  the  conservation  of  species 
equation.  The  equation  set  for  the  solid  phase  also  included 
source  terms,  describing  the  production  of  electrons  and  the 
generation  of  heat  due  to  the  electro-chemical  and  hetero¬ 
geneous  reactions.  Heat  can  also  be  transferred  between  the 
fluid  and  solid  phases,  through  a  source  term  representing 
convective  heat  transfer. 

The  equation  sets  for  the  fluid  and  solid  phases  are  in  a 
conservative  form.  Therefore,  future  work  will  involve  ap¬ 
plying  methods  from  computational  fluid  dynamics  (CFD) 
to  solve  the  equation  sets. 
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